!**********************************************************************
!*     F(I,J)  - начальный массив функции тока
!*     F1(I,J) - конечный массив функции тока
!*     Q1(I,J) - массив вихря
!*     S  - фиктивный шаг по времени
!*     EM - максимальная погрешность
!*     EO - относительная погрешность
!**********************************************************************
      SUBROUTINE PVR ( F, F1, Q1 )
      IMPLICIT REAL*8 ( A - Z )
      INTEGER I, J, L, N1, N2
      DIMENSION F(N1,N2), F1(N1,N2), Q1(N1,N2)
      COMMON /BLK/ S, HX, HY, N1, N2, EP
!******   Вычисление констант 
	!print *, "S=", S
      A1 = S/(HX*HX)
      A2 = S/(HY*HY)
!******   Начало итераций
      DO L = 1, 32000
        EM = 0.0d0
        
C$OMP PARALLEL
C$OMP DO PRIVATE(J) SHARED (F1, F, A1, A2, S) REDUCTION(MAX:EM)
        DO I = 2, N1 - 1
          DO J = 2, N2 - 1
            F1(I,J) = F(I,J) + A1*( F(I+1,J) - 2.0*F(I,J) + F(I-1,J) ) +
     &         A2*( F(I,J+1) - 2.0*F(I,J) + F(I,J-1) ) - S*Q1(I,J)
             EM = MAX( ABS( F1(I,J) - F(I,J) ), EM )
          ENDDO
        ENDDO
C$OMP END DO   
C$OMP END PARALLEL
      
      FS = 0.0d0
      
      DO I = 2, N1 - 1      
        DO J = 2, N2 - 1 
          FS = ABS( F1(I,J) ) + FS
        ENDDO
      ENDDO
      
      IF ( FS.LE.0.0d0 ) THEN
        PRINT *, ' Method PVR1 doesn^t work '
        STOP
      ENDIF
      
      EO = EM*DBLE( (N1-2)*(N2-2) )/FS
      !PRINT 120, F(4,9) 
120   FORMAT ( ' F(4,9) = ', 1P1E13.6 )

      IF ( EP.GT.EO ) EXIT
      
      DO I = 2, N1 - 1
        DO J = 2, N2 - 1 
            F(I,J) = F1(I,J)
        ENDDO
      ENDDO
      
      ENDDO
      
      RETURN
      END
